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Abstract 

We  develop  and  demonstrate  the  capability  of  a  high  order  accurate  finite  dif¬ 
ference  weighted  essentially  non-oscillatory  (WENO)  solver  for  the  direct  numerical 
simulation  of  transients  for  a  two  space  dimensional  Boltzmann  transport  equation 
(BTE)  coupled  with  the  Poisson  equation  modeling  semiconductor  devices  such  as  the 
MESFET  and  MOSFET.  We  compare  the  simulation  results  with  those  obtained  by  a 
direct  simulation  Monte  Carlo  (DSMC)  solver  for  the  same  geometry.  The  main  goal 
of  this  work  is  to  benchmark  and  clarify  the  implementation  of  boundary  conditions 
for  both,  deterministic  and  Monte  Carlo  numerical  schemes  modeling  these  devices, 
to  explain  the  boundary  singularities  for  both  the  electric  field  and  mean  velocities 
associated  to  the  solution  of  the  transport  equation,  and  to  demonstrate  the  over¬ 
all  excellent  behavior  of  the  deterministic  code  through  the  good  agreement  between 
the  Monte  Carlo  results  and  the  coarse  grid  results  of  the  deterministic  WENO-BTE 
scheme. 


Keywords:  Weighted  Essentially  Non-Oscillatory  (WENO)  schemes;  Boltzmann  Tran¬ 
sport  Equation  (BTE);  semiconductor  device  simulation;  MESFET;  MOSFET;  Direct  Sim¬ 
ulation  Monte  Carlo  (DSMC). 
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1  Introduction 


Electron  transport  in  semiconductors,  with  space  inhomogeneities  at  a  tenth  micron  scale 
orders  under  long  range  interactions  due  to  strong  applied  bias,  needs  to  be  modelled  by 
meso-scale  statistical  models.  One  such  type  of  models  is  the  semi-classical  approximation 
given  by  the  Boltzmann  Transport  Equation  (BTE)  for  semiconductors: 

^  +  ivk£  •  V,/  -  •  Vk/  =  Q(f).  (1.1) 

Following  standard  notation  for  statistical  descriptions  of  electron  transport  [21,  25],  / 
represents  the  electron  probability  density  function  (pdf)  in  phase  space  k  at  the  physical 
location  x  and  time  t.  Physical  constants  h  and  e  are  the  Planck  constant  divided  by  2ir 
and  the  positive  electric  charge,  respectively.  The  energy-band  function  £  is  given  by  a 
non-negative  continuous  function  of  the  form 


in  the  Kane  non-parabolic  band  model.  Here,  m*  is  the  effective  mass  and  a  is  the  non 
parabolicity  factor.  Thus,  by  setting  a  —  0  in  equation  (1.2)  the  model  is  reduced  to  the 
widely  used  parabolic  approximation. 

Our  model  takes  into  account  the  electrostatics  produced  by  the  electrons  and  the 
dopants  in  the  semiconductor.  Therefore,  the  electric  held  E  is  self  consistently  computed 
by 

Vx  [er(x)  VXK]  =  —  [p(£,x)  -  N D (x)] ,  (1.3) 

Co 

E  =  — VXK,  (1.4) 

where  eo  is  the  dielectric  constant  in  a  vacuum,  er(x)  labels  the  relative  dielectric  function 
depending  on  the  material,  p(f,x)  =  x,  k)  dk  is  the  electron  density,  iYo(x)  is  the 

doping  and  V  is  the  electric  potential.  Eqs.  (1.1)-(1.3)-(1.4)  give  the  Boltzmann-Poisson 
system  for  electron  transport  in  semiconductors. 

Regarding  the  collision  term  Q(f),  the  most  important  interactions  in  Si  are  due  to 
scattering  with  lattice  vibrations  of  the  crystal,  and  particularly,  those  due  to  the  acoustic 
and  optical  non-polar  modes.  Thus,  the  collision  term  will  read  as 

Q(f)(t,x,k)  =  f  [S(k',k)/(f,x,k') -S(k,k')/(£,x,k)]  dk',  (1.5) 

Jm3 
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where  the  kernel  S  is  defined  by 

S(  k,k')  =  X0(k,k')5(£(k')-£(k)) 

+  K (k,  k')  [(n?  +  l)5(e(k/)  -  e(k)  +  ftwp)  +  n(?5(e(k/)  -  e(k)  -  Hlop)\  .(1.6) 

Here,  we  have  taken  into  account  acoustic  phonon  scattering  in  the  clastic  approximation 
and  optical  non-polar  phonon  scattering  with  a  single  frequency  Hujp.  The  phonon  occupa¬ 
tion  factor  ng  is  given  by 

= [exp  (M)  -  *]  • 

where  ks  is  the  Boltzmann  constant  and  Tl  is  the  constant  lattice  temperature.  The  symbol 
5  indicates  the  usual  Dirac  distribution. 

From  the  distribution  function  /  one  evaluates  its  moments,  namely,  density  p(i,  x), 


mean  velocity 

1 

[  T(Vke)/(*,x,k)dk, 

'R3  n 

u(f,  x)  = 

1 

(1.7) 

p(f,x)  J 

and  mean  energy 

-1 

£(f,x) 

1 

-  /  e(k)/(t,x,k)  dk. 

)  Jr3 

(1.8) 

In  this  paper,  we  continue  our  research  to  find  suitable  deterministic  solvers  for  the 
Boltzmann-Poisson  system  of  semiconductors.  In  particular,  we  demonstrate  the  capabi¬ 
lity  of  our  WENO-BTE  scheme  [3]  in  giving  reliable  results  for  2D  realistic  devices.  We 
consider  MESFET  and  MOSFET  devices  which  have  served  as  benchmarks  for  several 
other  approximating  models  to  device  transport  simulation,  such  as  hydrodynamic,  drift- 
diffusion,  and  energy  transport  solvers  that  have  typically  been  compared  to  Monte  Carlo 
results  for  solving  the  BTE  [1,  7].  We  refer  to  [8,  6,  9]  for  previous  and  alternative  attempts 
for  developing  ID  deterministic  solvers  for  BTE. 

We  show  that  the  results  of  our  deterministic  solver  WENO-BTE  agrees  remarkably  well 
with  the  ones  obtained  by  the  Monte  Carlo  approach  in  the  areas  of  the  device  in  which 
there  are  enough  particles  to  obtain  sufficient  statistics  for  particles.  The  main  advantage 
of  the  WENO-BTE  solver  comes  from  its  suitability  to  describe  almost-empty  areas  of  the 
device,  where  it  is  quite  hard  for  Monte  Carlo  techniques  to  obtain  relevant  information  on 
density,  mean  velocity  and  energy.  These  almost-empty  regions  are  typically  caused  by  the 
strong  helds  repelling  the  electrons  near  the  gate/drain  regions. 

In  addition,  deterministic  solvers  provide  a  tool  to  compute  accurate  transient  calcula¬ 
tions  of  strong  non-statistical  equilibrium  states  as  was  already  shown  in  [3]  for  one  space 
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dimensional  solvers  simulating  diodes.  Our  simulations  for  the  two  dimensional  case  indicate 
that  one  observes  stabilization  of  moments  in  time  and  thus,  we  use  small  time  variations  of 
the  density,  current  and  energy  as  the  code  stopping  criteria.  A  parallelization  of  the  code 
for  a  fine  grid  has  produced  numerical  evidence  ensuring  that  stable  states  are  obtained  in 
higher  dimensions  as  well  [20]  (see  the  2D  MESFET  transient  simulations  by  parallelized 
WENO-Boltzmann  schemes  in  http :  / /kinetic  .mat  .uab .  es/~carrillo/galleryl  .html). 


In  addition,  deterministic  solvers  allow  for  the  explicit  computation  of  the  probabil¬ 
ity  density  function  solution  /  to  the  BTE  transport  model,  and  thus,  all  moments  are 
computed  with  a  better  accuracy  compared  to  Monte  Carlo  codes.  Therefore,  closure  as¬ 
sumptions  typically  used  in  hydrodynamic  systems  can  be  validated  numerically  with  our 
tools,  as  well  as  reliable  code  hybridization  approaches  linking  kinetic  to  hydrodynamics 
under  long  range  interactions,  by  identifying  appropriate  stiff  regimes  where  hydrodynamic 
asymptotics  are  achieved.  Consequently,  our  deterministic  solvers  can  be  considered  as 
the  real  benchmark  for  moment-based  transport  solvers:  hydrodynamic  [1],  drift-diffusion, 
energy-transport  models  [7]  and  so  on;  instead  of  Monte  Carlo  results. 

Following  the  main  ideas  of  [19,  2,  3],  we  have  introduced  in  our  previous  work  [4] 
a  suitable  dimensionless  formulation  by  a  change  of  variables  in  the  wave  vector  k  that 
translates  the  original  Boltzmann  equation  (1.1)  in  2D  into  the  following  system:  we  use 
the  coordinate  transformation  for  the  momentum  vector  k 


k  = 


\Z2m*kBTL 

h 


+  cxkuj  i^y,  \/l  —  y2  cos  4>,  \/l  —  y2  sin  0  j  (1.9) 


where  olk  =  kBTLa,  u>  >  0  is  a  dimensionless  energy,  y  G  [—1, 1]  is  the  cosine  of  the  angle 
with  respect  to  the  x-axis  and  0  €  [0,  27r]  is  the  azimuthal  angle,  respectively.  We  consider 
a  two  dimensional  problem  in  the  physical  space  and  we  define 


$(f,  x ,  y,  ix,  y,  0)  =  s(u)f(t,  x,  y,  u,  y,  0)  , 

where  s{u)  =  \/ix(  1  +  ~a^uj)(l  +  2aKu)  is,  apart  a  dimensional  constant  factor,  the  Jacobian 
of  the  transformation  of  coordinates  (1.9)  and  /  corresponds  to  the  distribution  function 
/  after  the  change  of  variable  (1.9).  The  use  of  the  dimensionless  unknown  $  instead  of 
the  distribution  function  /  allows  us  to  write  the  dimensionless  Boltzmann  equation  in 
conservative  form  as 

W  +  ^(“1*)  +  +  Ia-4’)  +  =  sMC(4>).  (1.10) 
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where  the  new  collision  operator  C,(<h)  and  the  flux  functions  a,  were  computed  in  [4]  and 
are  repeated  here  for  the  sake  of  the  reader:  the  flux  functions  a*  are  given  by 


ai(ca)  = 


fis(ui ) 


a2(  u,n,4>)  = 


U  (1  +  2 aKcu)2 

1  \Jl  —  n2s(u)  cos  0 


i  %  i  y  ^  1 4>) 

O' 4(^7  V7  ^ 7  AL  0) 

a5(t,x,y,uj,n,(f))  = 


U  (1  +  2aKu;)2 

1  2s(u) 


U  (1  +  2aKa;)2 
1  1  H~  2aKco 


VT 


U  s(u>) 

Ey(t ,  x,  y )  sin  0  1  +  2aKu; 

t*  y/1  ~  /l2  S{(J) 


Ex(t,  x,  y)n  +  Ey(t,x,y)^l  -  y2  cos  0 

Ex(t,x,y)y/l  -  y2  -  Ey(t,x,y)n cos0 


and  the  dimensionless  collision  operator  C(<h)  is  given  by 
C($)(t,x,y,uj,y,(j))  = 

27tT  J0  /  Xl  V ’ ^  ^  +  x,y,LJ  +  a>  ll>1  ^  a>  /A  0')]#/(V 


1 


s(u;)  t 


■[/3s (a;)  +  as(cu  —  a)  +  s(cu  +  a)]<h(t,  x,  y,  u,  y,  0) . 


The  constant  parameters  f*,  a  and  (3  depend  on  the  scattering  mechanisms.  The  constant 
a  is  the  dimensionless  phonon  energy  and  the  energy  mesh  has  to  be  chosen  properly  for 
computational  simplification  of  the  collision  kernel  computation.  For  all  details  and  physical 
constants  we  refer  to  works  [2,  3,  4], 


Let  us  remark  that  preliminary  results  of  2D  Si  MESFET  devices  were  presented  in  [4] 
and  that  very  fine  grid  results  of  a  2D  Si  MESFET  and  its  transients  were  reported  in  [20] 
based  on  a  parallelization  of  the  numerical  code  we  will  present  in  this  work.  Our  main 
objective  in  this  work  is  to  clarify  the  implementation  of  boundary  conditions  for  the  deter¬ 
ministic  and  for  the  Monte  Carlo  numerical  scheme,  explain  the  appearance  of  singularities 
in  the  electric  field  and  demonstrate  the  overall  excellent  behavior  of  the  deterministic  code 
through  the  good  agreement  between  the  Monte  Carlo  results  and  the  coarse  grid  results 
of  our  deterministic  WENO-BTE  scheme. 


The  structure  of  the  paper  is  as  follows.  In  the  next  section,  we  describe  the  numerical 
implementation  of  the  kinetic  boundary  conditions  used  within  the  deterministic  compu¬ 
tation  and  we  summarize  the  main  characteristics  of  the  deterministic  scheme  previously 
discussed  in  [4],  We  also  discuss  the  deterministic  simulation  of  the  2D  MESFET  and 
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MOSFET,  including  its  efficiency,  accuracy  and  performance  in  this  section.  In  Section  3, 
we  give  a  presentation  of  the  Monte  Carlo  method  that  has  been  used  to  compare  to  the 
deterministic  scheme.  Finally,  we  draw  our  conclusions  in  Section  4.  Our  presentation  in 
Sections  2  and  3  will  be  based  on  the  particular  configuration  of  the  MESFET  device  [16], 
since  it  contains  all  the  main  difficulties  we  will  deal  with,  and  we  will  present  the  results 
for  the  MOSFET  [25]  as  a  mere  application  of  all  the  knowledge  and  difficulties  faced  and 
solved  in  the  MESFET  case. 


2  2D  deterministic  device  simulations 

We  describe  the  electron  transport  in  a  MESFET  by  solving  the  initial  boundary  value 
problem  associated  to  equations  (1.3),  (1.4)  and  (1.10)  in  a  spatial  rectangular  domain 
denoted  by  R  =  [0,  x0]  x  [0,2/o]  (see  Figure  1)  and  with  all  R3  as  the  phase-velocity  space. 
Our  first  objective  in  this  section  is  to  describe  in  detail  the  boundary  conditions  and  its 
numerical  implementation. 

2.1  Boundary  conditions  for  the  BTE-Poisson  system 

We  describe  now  the  boundary  conditions  we  have  used  to  simulate  both  MESFET  and 
MOSFET  geometries  by  following  classical  assumptions  for  the  modelling  of  two  dimen¬ 
sional  devices  [21,  25]. 

At  the  source  and  drain  contact  regions,  numerical  boundary  conditions  must  approxi¬ 
mate  neutral  charges. 

At  the  gate  area,  the  numerical  boundary  condition  yields  an  estimated  incoming  mean 
velocity  which  represents  a  transversal  electric  field  effect  clue  to  the  gate  contact  (see  the 
electric  held  in  Figure  3).  These  boundary  conditions  are  imposed  in  an  outside  buffer  layer 
above  the  contact  and  gate  areas.  As  in  classical  simulations  of  hyperbolic  systems  with 
methods  such  as  WENO  schemes  for  flux  reconstruction,  these  inflow  conditions  consist  of 
the  approximation  of  incoming  flow  on  characteristic  lines  that  cross  the  boundary  towards 
the  interior  of  the  domain.  That  means  these  conditions  are  consistent  with  the  well- 
posedness  of  the  initial  boundary  value  problem  of  the  corresponding  BTE-Poisson  system. 

Insulating  conditions  are  imposed  in  the  remainder  boundary  of  the  device. 
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Even  though  the  existence  of  solutions  to  the  initial-boundary  value  problem  related 
to  our  continuum  model  is  open,  our  simulation  indicates  that  the  solution  is  regular  in 
the  interior  of  the  domain  and  may  develop  boundary  singularities  at  the  junctions  of  the 
contacts  and  insulating  walls.  In  addition,  these  singularities  are  consistent  with  the  same 
type  of  behavior  for  boundary  problems  for  elliptic  and  parabolic  equations  with  boundary 
junctions  of  Dirichlet  and  Neumann  type  conditions,  which  are  described  below. 

As  for  the  boundary  conditions  associated  to  the  electrostatic  potential,  (i.e.  solutions 
of  the  Poisson  equation  (1.3)-(1.4)  ),  we  impose  classical  boundary  conditions  for  elliptic 
solvers  associated  to  two-dimensional  devices.  They  are  as  follows:  prescribed  potential 
(Dirichlet  type  conditions)  at  the  source,  drain,  and  gate  contacts,  and  insulating  (homo¬ 
geneous  Neumann  conditions)  on  the  remainder  of  the  boundary. 

Let  us  denote  by  <9R  =  <9Rd  U<9Rjv,  the  boundary  of  the  device  which  we  separate  into 
<9Rd  for  the  gates  and  contact  areas,  and  <9R n  for  the  insulated  areas.  In  particular  the 
top  boundary  (x,yo)]x  G  [0,xo]  contains  the  source,  drain  and  gate  contacts.  In  classical 
modelling  of  micro-electronic  devices,  the  doping  function  Ab)(x)  =  Nd[x,  y )  is  only  defined 
in  the  semi-conducting  domain  R.  It  is  assumed  to  be  a  piecewise  constant  function, 
with  jump  surfaces  separating  un-doped  regions,  denoted  by  n(x,  y),  where  the  background 
density  is  just  given  by  bulk  material  properties,  from  strongly  doped  regions,  denoted  by 
n+(x,  y),  where  density  of  background  electrons  has  been  increased  to  higher  concentrations 
than  the  ones  corresponding  to  bulk  material.  Source  and  drain  contact  boundary  regions 
must  satisfy  ND(x.)  =  n+(x,y),  that  is,  are  boundary  regions  located  on  areas  of  high 
background  concentration. 

More  specifically,  the  boundary  conditions  imposed  to  the  numerical  probability  density 
function  (pdf)  approximating  the  solution  of  the  BTE,  are  as  follows: 

Source  and  drain  contacts:  We  employ  a  buffer  layer  of  ghost  points  outside  both 
the  MESFET  and  MOSFET  domain  in  these  contact  areas.  For  {y  =  yo\  and  Ay  >  0, 

-  Nn(r) 

fit,  x,  y0  +  Ay,  k)  =  —7. — - - - f(t,  x,  y0,  k)  (2.1) 

/  f(t,x,y0,k!)dk! 

Jr 3 

where  N0(x)  =  ND(x,y0)  is  the  doping  profile  and  f(t,x,y0,  k)  denotes  the  evaluation  of 
the  numerical  solution  on  mesh  points  located  at  {t,  x,  y0,  k(cn,  y,  (f>)),  where  x  lies  in  those 
intervals  such  that  the  pair  ( x ,  yf)  is  either  a  point  in  the  source  or  drain  boundary  region. 

Linder  the  assumption  that  highly  doped  regions  n+(x,y)  are  slowly  varying  (i.e.  of 
very  small  total  variation  norm),  the  corresponding  small  Debye  length  asymptotics  yields 
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neutral  charges  away  from  the  endpoints  of  the  contact  boundary  regions.  Consequently, 
the  total  density  p(t,  x,  y0 )  (zero  order  moment)  and  its  variation  takes  asymptotically,  in 
Debye  length,  the  values  of  n+{x,y)  away  from  the  contact  endpoints.  This  is  the  well- 
known  limit  for  Drift  Diffusion  systems  [21],  and  we  observe  that  this  is  also  the  case  for 
the  kinetic  problem,  provided  the  kinetic  solution  satisfies  approximately  neutral  charge 
conditions  at  these  contact  boundaries. 

Thus,  under  the  assumption  that  n+(x,y )  is  sufficiently  large  and  that  \'Vyn+(x,y)\  is 
small,  the  Debye  asymptotics  near  the  source  and  drain  contacts  must  yield  p  ~  n+  and 
\Wyp\  ~  |V2/n+|,  so  p(x,y)  has  bounded  first  order  partial-?/  derivative  at  y0  for  (x,y0)  a 
point  in  the  interior  of  the  contact  boundary. 

In  particular,  a  simple  calculation  shows  that  condition  (2.1)  yields  neutral  charges  at 
these  source  and  drain  boundary  regions.  Indeed,  integrating  in  velocity  space  both  sides 
of  the  identity  (2.1),  for  p(t,x,y0 )  =  f(t,x,y0,k)  dk 

p(t,x,y0  +  Ay)  =  N0(x) 


or,  equivalently, 


pit,  x ,  yo )  «  N0{x,  ?/o )  -  Ay  Vyp(t,  x,  y0)  « 

N0(x,  2/0 )  -  Ay  VyN0(x,  y0 )  «  N0(x,  y0)  +  0(77)  (2.2) 

where  rj  =  rj(Ay,  A,  HiVollc1))  and  A  the  scaled  Debye  length  parameter  associated  to  the 
contact  area  for  highly  doped  regions. 


Gate  contact:  We  first  describe  the  boundary  condition  for  the  MESFET  simulation. 
Here  we  set  a  condition  that  yields  an  incoming  mean  velocity  entering  through  the  gate 
which  is  equivalent  to  decaying  density  in  space  and  time  due  to  the  conservation  of  charge. 
Thus  we  set 

fit ,  x,  2/0  +  Ay,  k)  =  nf(t,  x,  2/0,  k),  (2.3) 

then,  again  integrating  in  velocity  space  both  sides  of  the  identity  (2.3), 


\/yp{t,x,y0)  « 


K  —  1 

Ay 


p{t,x,y0). 


(2.4) 


In  fact,  boundary  condition  (2.4),  which  is  of  Robin  (or  mixed)  type,  is  the  natural 
simulation  of  a  gate  contact  in  an  oxide  region.  The  quantity  (k  —  1)/A y  must  reflect  the 
scale  of  the  electric  field  flux  across  the  interface  corresponding  to  the  gate  voltage.  This 


quantity  is  related  to  the  quotient  of  the  permittivity  and  thickness  of  the  gate  region  where 
charges  are  assumed  negligible  (oxide). 

In  order  to  explain  this  effect,  we  recall  classical  asymptotics  corresponding  to  thin  shell 
elliptic  and  parabolic  type  problems  with  Dirichlet  data  in  the  thin  outer  shell  domain  and 
a  core  larger-scale  domain,  where  classical  transmission  conditions  link  both  domains.  The 
two  dimensional  asymptotic  boundary  condition  for  the  Poisson  system  linking  the  oxide 
and  semiconductor  region  yielding  a  Robin  type  condition  was  rigorously  studied  in  [11,  12] 
and  is  given  by 

<25) 

where  r]  denotes  the  outer  unit-vector  to  the  interface  surface  given  by  the  intersection  of 
the  boundaries  to  the  oxide  and  the  semiconductor  regions.  The  value  Vgate  is  the  prescribed 
electrostatic  potential  at  the  gate  contact  and  eox  and  5  are  the  prescribed  oxide  permittivity 
and  thickness  respectively.  This  means  the  prescribed  gate  voltage  Vgate  is  on  the  edge  of 
the  eox-gate  at  5-distance  from  the  semiconducting  domain. 

Therefore,  in  the  asymptotic  limit  where  0  <  eox/5  =  C  <  oo  as  both  eox  — >  0  and 
5  — >  0,  condition  (2.5)  becomes 

— E  ■  jj  —  W  ■  7/  —  h—  C-f(V  -  Vgate).  (2.6) 

which  is  a  Robin  type  boundary  condition  for  the  Poisson  equation  at  the  oxide-semiconductor 
interface.  This  asymptotic  boundary  condition,  formally  derived  for  drift-diffusion  equa¬ 
tions  by  the  computational  electronic  community,  have  been  widely  used  in  device  simula¬ 
tions  (see  [22]  and  references  therein). 

Hence,  in  our  present  setting,  one  can  interpret  the  relation  (2.4)  as  a  choice  of  an 
electric  held  flux  rate  across  the  gate  interface,  where  Ay  is  chosen  a  multiple  of  the 
oxide  thickness,  Ay  &  k  5  and  then  n  is  chosen 

1  +  —  R  te),  (2.7) 

27t  e 

when  V  is  available. 

Then  condition  (2.4)  becomes,  after  combining  with  (2.6)  and  (2.7), 

K  —  1 

Vyp(t,x,y0)  «  -Eyp(t,x,y0) ,  ~Ey  =  -E  •  rj  =  (2.8) 

which  coincides  with  a  classical  macroscopic  boundary  condition  for  a  density  p  that  evolves 
according  to  a  drift- diffusion  equation.  In  general,  the  value  Vgate  is  a  prescribed  gate  voltage 
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and  V  is  the  computed  value  of  the  voltage  in  the  gate  interface  for  the  Poisson  solver  in  a 
semiconductor-oxide  region  with  boundary  conditions. 

We  remark  that  our  MESFET  simulations  do  not  solve  Poisson  equation  (1.3)  using  the 
asymptotic  boundary  condition  (2.6)  for  small  (5-thick  gate/oxide  region.  The  main  reason 
is  that  we  are  neglecting  the  oxide  region,  and  consequently  we  need  to  impose  the  value  of 
the  potential  at  the  gate  to  be  Vgaie  and  fix  some  boundary  condition  for  the  distribution 
function  /.  Instead,  in  order  to  compensate  the  held  flux  across  the  gate  interface,  we 
having  fixed  k  =  0.5  and  Ay  a  unit  of  mesh  size.  This  choice  is  asymptotically  equivalent 
to  have  fixed  V  on  the  MESFET  gate  region  depending  on  Vgate  and  the  oxide  thickness  6 
in  such  a  way  that  Ay  =  R  5,  so  that 

7"  —  k  1  ~  eA^(V  -  Vgate)  or  equivalently  V  «  Egate - ^-4r-  (2.9) 

Z  Zll  6  6 qX  ft 

In  particular,  the  choice  of  the  parameter  k  yields  an  approximation  that  may  modify 
the  effective  held  hux  across  the  gate  contact.  In  fact  this  choice  may  be  interpreted  as  a 
calibrating  parameter  to  evaluate  the  strength  of  the  electric  held  hux  across  the  gate/oxide 
interface  when  Dirichlet  data  V^ate  has  been  prescribed  instead  of  either  the  Robin  type  data 
(2.6)  or  the  more  accurate  simulation  of  the  Poisson  solver  in  the  whole  oxide/semiconductor 
domain. 

Finally,  for  the  MOSFET  simulation  the  gate  contact,  we  recall  that  the  Poisson  solver 
imposes  a  Dirichlet  condition  by  prescribing  Vgate  and  since  the  gate  contact  is  outside 
the  domain  of  the  transport  equation,  no  condition  is  needed  for  the  computation  of  the 
probability  density  function  pdf. 

Insulating  walls:  In  the  remainder  of  the  boundary  region,  we  impose  classical  clastic 
specular  boundary  reflection.  For  x  e  <9Rjv 

f(t,  x,  k)  =  f(t,  x,  k*)  for  k  •  ?7(x)  <  0  with  k*  =  k  —  2k  •  77(x)  (2-10) 

dV 

-(,x)=0.  (2. 

Oxide/semiconductor  MOSFET  interface:  The  Poisson  equation  is  posed  with 
neutral  charges  in  the  oxide  region  and  (1.3)  in  the  semiconducting  region.  Condition  (2.10) 
is  imposed  on  the  oxide/semiconductor  interface  and  the  electric  held  on  this  interface  is 
computed  from  the  Poisson  solver. 
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2.2  Singular  boundary  behavior. 

2.2.1  MESFET  geometry. 


At  the  boundary  region  y0  =  0.2 gm,  we  have  contact/insulating  and  insulating/gate 
changes  of  the  boundary  condition.  It  is  well-known  that  solutions  of  the  Poisson  equation 
develop  singularities  under  these  conditions.  Namely,  when  solving  the  Poisson  equation 
in  a  bounded  domain,  with  a  Lipschitz  boundary  and  a  boundary  condition  that  changes 
type  from  Dirichlet  to  homogeneous  Neumann  or  Robin  type  data  on  both  sides  of  a  conical 
wedge,  with  an  angle  0  <  d  <  2i r,  whose  vertex  is  at  the  point  of  data  discontinuity,  the 
solution  will  develop  a  singularity  at  the  vertex  point  depending  on  the  angle  R. 

More  specifically,  for  an  elliptic  problem  in  a  domain  with  data-type  change,  the  solution 
remains  bounded  with  a  radial  behavior  of  a  power  less  than  one  if  vr/2  <  R  <  2n.  Its 
gradient  then  becomes  singular  (see  for  instance  Grisvard  [14]  for  a  survey  on  boundary 
value  problems  of  elliptic  PDE’s  in  a  non-convex  domain). 

In  this  case,  it  is  possible  to  prove  the  corresponding  asymptotic  behavior  is  controlled 
from  above  and  below  by  functions  with  bounded  Laplacian  whose  radial  component  is  of 
order  0(W2l?).  More  precisely,  a  local  control  function  (barrier)  to  solutions  of  Laplace’s 
equations  with  bounded  right  hand  side  can  be  constructed  by  taking  functions  of  the  form 
rT  sin(r6)  +  CQ  with  r  = 

In  the  particular  case  of  the  jump  of  a  boundary  condition  in  a  Rat  boundary  surface 
(zero  curvature),  such  as  in  the  boundary  given  by  y$  =  0.2,  the  angle  is  L?  =  7 r  which  makes 
the  electric  potential  V  =  r1//2  sin(6/2)  +  C0.  Then  its  gradients  (i.e.  the  components  of 
electric  field)  develop  a  singularity  of  the  order  W  ~  Mr-1/2.  Similar  local  asymptotic 
behavior  occurs  for  parabolic  systems  in  bounded  domains  with  discontinuous  boundary 
data  operators. 

We  remark  that  our  calculations  resolve  this  singular  boundary  behavior  as  is  de- 
mostrated  in  Figures  2  and  3. 

Consequently,  elliptic  and  parabolic  descriptions  of  the  density  flow  at  the  drift-diffusion 
level  remain  with  the  same  singular  boundary  behavior,  even  though  mobilities  functions  are 
saturated  quantities;  that  is,  they  are  bounded  functions  of  the  magnitude  of  the  electric 
field  [12].  In  fact,  these  conditions  have  been  used  for  two  dimensional  drift- diffusion 
simulations  of  MESFET  devices  and  this  singular  boundary  behavior  was  obtained  [5]. 
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In  view  of  this  boundary  asymptotic  behavior  at  the  macroscopic  level  and  the  numerical 
results  we  obtain  for  the  deterministic  solutions  of  the  BTE-Poisson  system  and  its  moments 
(see  Figures  2  and  3),  we  conjecture  that  the  moments  of  the  particle  distribution  function 
/(t,x,  k)  have  the  same  spatial  asymptotic  behavior  as  the  density  solution  of  the  drift 
diffusion  problem  at  the  boundary  points  with  a  discontinuity  in  the  data.  In  other  words, 
the  collision  mechanism  preserves  the  spatial  regularity  of  the  pdf  as  though  its  zeroth 
moment  would  satisfy  the  drift- diffusion  equations,  even  though  the  numerical  evidence 
indicates  that  the  average  of  the  pdf  (i.e.  density)  will  not  evolve  according  to  such  a 
simple  macroscopic  model. 

Based  in  our  calculations  in  Figures  2  and  7,  we  conjecture  the  density  p(f,x)  develops  a 
contact  local  singularity  of  the  same  asymptotic  order  as  the  electrostatic  potential  V  (t,  x) 
at  (t,x,yo).  In  particular,  under  the  assumption  that  the  first  and  second  moments  of  the 
/(t,x,  k)  are  bounded  in  (f,x),  local  behavior  of  the  mean  velocity  w(t,x)  defined  in  (1.7) 
compares  to  the  one  of  the  electric  field  E(t,  x),  and  thus 

\u(t,x,y0)\  ~  Mr^1/2  (2-12) 

with  r  the  distance  from  (x,y)  to  the  vertex  where  the  data  jumps.  Similarly,  the  mean 
energy  S(t,x,y0),  defined  in  (1.8),  develops  a  singularity  of  the  same  order  as  the  mean 
velocity  in  (2.12). 

2.2.2  MOSFET  geometry. 

In  this  particular  case,  the  electron  pdf  transport  is  modelled  only  in  the  semiconducting 
region  with  insulation  conditions  along  the  oxide/semiconductor  interface.  The  first  five 
moments  are  expected  to  be  bounded  along  the  interface  (see  Figure  11).  The  only  points 
where  singularities  generating  singular  electric  fields  could  arise  for  the  Poisson  solver  are 
the  points  where  the  n+  —n  junctions  meet  the  semiconductor-oxide  interface.  We  observe 
that,  in  fact,  the  maximum  value  for  the  energy  (second  moment  for  the  pdf)  occurs  at  the 
n  —  n+  drain  junction. 

Indeed,  it  has  been  shown  in  [10]  that  all  second  derivatives  of  the  potential  will  remain 
bounded  at  those  intersections,  provided  the  boundary  condition  along  the  interface  are 
of  transmission  or  insulating  type  for  a  drift-diffusion  model.  From  our  calculations  (see 
Figure  11),  we  expect  similar  behavior  for  the  spatial  regularity  of  the  moments  associated 
to  the  kinetic  solution  as  well,  which  in  these  case,  would  yield  bounded  mean  velocities 
implying  bounded  energy  throughout  the  semiconducting  device  domain. 
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2.3  The  WENO-BTE  solver 


We  adopt  the  same  fifth  order  finite  difference  WENO  scheme  coupled  with  a  third  order 
TVD  Runge-Kutta  time  discretization  to  solve  (1.10)  as  that  used  in  our  previous  work 
[3].  The  fifth  order  finite  difference  WENO  scheme  was  first  designed  in  [17],  and  the  third 
order  TVD  Runge-Kutta  time  discretization  was  first  proposed  in  [24],  In  our  calculation, 
the  computational  domain  in  each  of  the  Eve  directions  x,  y,  to,  fi  and  0  is  discretized  into 
a  uniform  grid,  with  1  +  Nx,  1  +  Ny,  Nu,  N y  and  points  in  these  directions  respectively. 
Our  numerical  method  also  allows  the  use  of  non-uniform  but  smooth  grids  in  order  to 
concentrate  points  near  regions  where  the  solution  has  larger  variations.  In  the  x  and  y 
directions,  the  computational  boundaries  are  located  at  grid  points.  In  the  other  three 
directions,  the  computational  boundaries  are  located  in  the  middle  of  two  grid  points  in 
order  to  efficiently  implement  symmetry  type  boundary  conditions  to  maintain  conservation. 

The  approximations  to  the  point  values  of  the  solution  &(tn,Xi,yj,Uk,  ni,(j>m),  denoted 
by  •  k  i  m,  are  obtained  with  a  dimension  by  dimension  (not  dimension  splitting)  approxi¬ 
mation  to  the  spatial  derivatives  using  the  fifth  order  WENO  method  in  [17].  For  example, 
when  approximating  ^  (ai<h),  the  other  four  variables  y,  to,  /x  and  0  are  fixed  and  the 
approximation  is  performed  along  the  x  line: 

d  1 

(Ol(^fc)  ,  Xj,  yj ,  UJk,  fl[,  0m))  ~  ^ ~  (j^i+1/2  hi— 1/2^ 

where  the  numerical  flux  hi+ 1/2  is  obtained  with  a  one  dimensional  procedure.  We  will  not 
repeat  the  details  of  the  computation  for  this  numerical  flux  here  and  refer  the  reader  to 
[3,  17,  23], 

Next,  we  briefly  describe  the  numerical  implementation  of  the  boundary  conditions.  See 
the  previous  subsection  for  more  details  of  some  of  these  boundary  conditions. 

•  At  x  —  0  and  x  =  Xq,  we  implement  a  reflective  boundary  condition  by  prescribing 
the  values  of  the  ghost  points  using  for  x  >  0, 

$(t,  -x,  y,  to,  n,  0)  =  $(t,  x,  y,  to,  -/1,  0), 

<h(f,  Xq  +  x,  y,  to,  /x,  0)  =  $(t,  x0  -  x,  y,  to,  -/I,  0). 

At  y  =  0,  we  implement  a  similar  reflective  boundary  condition  using  for  y  >  0, 

$(f,  x,  -y,  to,  fi,  0)  =  $(t,  x,  y,  to,  /i,  2t t  —  0), 

for  the  ghost  points. 
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•  For  the  top  boundary  y  =  y0,  the  boundary  conditions  are  explained  in  the  previous 
subsection. 

•  At  u  =  0,  the  boundary  condition  is  given  to  take  into  consideration  that  this  is  not 
really  a  physical  boundary;  however,  a  “ghost  point”  at  (o J,y,4>)  for  a  negative  uj  is 
actually  a  physical  point  at  (—a;,  —y,  2n  —  0).  We  thus  take  the  boundary  condition 
for  uj  >  0  to  be, 


$(t,  x,  y ,  -uj,  y,  0)  =  $(t,  x,  y,  uj,  -y,  2i t  -  </>). 

An  averaging  is  then  performed  for  the  fluxes  for  uj  =  0  to  enforce  conservation. 

•  The  boundary  region  uj  =  ujmax  is  taken  as  an  outflow  boundary  with  a  Neumann 
type  condition 


x,  y ,  uj ,  y,  0)  &(t,  x,  y,  cumax,  y,  0) ,  uj  j>  cumax 

and  the  explicit  enforcement  of  /pv^+i/2  =  0  for  the  last  numerical  flux  in  order  to 
enforce  conservation  of  mass. 

•  At  y  —  —  1  and  y  —  1,  we  take  the  boundary  conditions  for  y  >  0, 

$(t,  x,  y,  uj,  -1  -  y,  0)  =  $(t,  x,  y,  uj, -1  +  y,</>), 

$(t,  x,  y,  w,  1  +  y,  0)  =  $(i,  x,  y,  uj,  1  -  y,  0), 

motivated  by  the  physical  meaning  of  y  as  the  cosine  of  the  angle  to  the  2- axis.  We 
also  explicitly  impose  h0  =  fiN^+1/2  =  0  for  the  first  and  last  numerical  fluxes  in  order 
to  enforce  conservation  of  mass. 

•  At  0  =  0  and  0  =  2n,  we  take  the  boundary  conditions  for  0  >  0, 

$(t,  x,  y,  uj,  y,  -0)  =  $(t,  x,  y,  uj,  y,  0), 
x,  y,  uj,  y,  2vr  +  0)  =  <f>(t,  x,  y,  uj,  y,  2vr  -  0), 

motivated  by  the  physical  meaning  of  0  as  the  angle  around  the  2-axis.  We  also 
explicitly  impose  h0  =  Av0+ 1/2  =  0  for  the  first  and  last  numerical  fluxes  in  order  to 
enforce  conservation  of  mass. 
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An  integration  to  compute  the  various  moments  using  the  pdf  is  approximated  by  the 
composite  mid-point  rule: 

/*+oo  /*1  p2tt 

/  du)  I  dfi  I  d(j)  b(u>,  y,(j))$(t,x,y,u>,  y,(j)) 

Jo  J- 1  Jo 

N<t>  Nu. 

-  EEE  /h>  0m)^=*(C  %i  Vi  ^ki  d'li  0m)  AcU A/l A0 

m=  1  i=l  fc=l 

which  is  at  least  second  order  accurate.  This  condition  is  even  more  accurate  when  <f>  and 
some  of  its  derivatives  are  zero  at  the  boundaries,  and,  in  fact,  it  is  infinite  order  accurate  or 
spectrally  accurate  for  C°°  functions  with  compact  support.  Notice  that  mass  conservation 
is  enforced  with  this  integration  formula  and  boundary  conditions  indicated  previously. 


The  Poisson  equation  (1.3)  for  the  potential  V  and  the  electric  field  (1.4)  are  solved 
by  the  standard  central  difference  scheme,  with  the  given  V0as  boundary  conditions  at  the 
source,  drain  and  gate. 

The  time  discretization  is  by  the  third  order  TYD  Runge-Kutta  method  [24], 


Finally,  due  to  the  explicit  character  of  the  time  evolution  solver  we  need  to  impose  the 
standard  CFL  condition  for  these  RK  schemes.  That  is, 


At  <  CFL 


Ax 


+ 


Ay 


+ 


Au 


+ 


Ay 


+ 


A  0 


max(|ai|)  max(|a2|)  max(|a3|)  max(|a4|)  max(|a5|) 


(2.13) 


where  the  maximum  is  taken  over  the  grid  points  and  we  use  CFL  =  0.6  in  the  actual 
calculation.  Note  that  we  are  taking  max(|a4|)  over  a  grid  without  u>  =  0.  For  the  cases 
we  have  run,  the  source  collision  term  is  not  stiff  enough  to  change  the  CFL  condition 
determined  by  the  approximations  to  the  spatial  derivatives.  This  CFL  condition  at  first 
glance  may  seem  severe,  however  some  practical  cases  shown  in  next  sections  demonstrate 
this  is  not  the  case  for  coarse  grids.  This  CFL  condition  does  become  more  severe  due  to 
the  appearance  of  the  singularities  in  the  force  held  for  more  refined  meshes.  Therefore, 
although  the  CPLT  time  for  coarse  grids  is  comparable  with  DSMC  simulations,  fine  grid 
high  accurate  deterministic  results  are  more  costly. 


2.4  Deterministic  simulation  of  a  2D  Si  MESFET 

In  Figures  2  and  3  we  show  the  results  of  the  macroscopic  quantities  for  the  MESFET  when 
a  grid  of  43  x  25  points  is  used  in  space,  80  points  in  the  energy  u,  and  12  points  in  each 
of  the  angular  y  and  0  directions  (grid  A).  The  doping  is  3  x  101'  cm-3  in  the  n+  zone 
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and  1017  cm-3  in  the  n  zone.  We  assume  applied  biases  of  Vsource  =  OW,  Vgate  =  —0.81/ 
and  Vdrain  =  1.1/.  The  velocity  and  electric  vector  fields  give  a  clear  picture  of  the  electron 
transport  in  the  device. 

We  demonstrate  numerically  that  singularities  appear  in  the  electric  field  and  the  square- 
root  behavior  of  the  electric  potential  is  clearly  defined.  The  same  behavior  seems  to  be 
the  case  for  the  density  and  the  velocity  held  but  the  spatial  resolution  is  not  good  enough 
for  a  complete  numerical  answer  to  this  point.  The  appearance  of  these  singularities  is 
a  difficulty  for  the  numerical  resolution  of  this  problem,  since  the  CFL  condition  (2.13) 
becomes  more  restrictive  due  to  the  singularity  of  the  electric  held  as  the  grid  gets  finer. 

We  observe  highly  accurate  non-oscillatory  behavior  near  junctions  and  expected  high 
channel  anisotropic  energy  near  the  gate  drain  region,  where  the  two-dimensional  effects  are 
most  signiheant.  These  assertions  are  also  based  on  the  results  of  the  probability  density 
functions  in  Figure  4.  These  pdf’s  have  been  obtained  by  averaging  the  computed  values 
of  <f>  over  the  angle  <fi.  We  observe  the  overall  smoothness  of  the  pdf  <f>  function  even  near 
the  gate-insulated  data-type  change  regions. 

We  also  show  in  Figures  5  and  6  the  discrepancies  between  the  results  obtained  with 
grid  A  and  the  results  of  two  other  grids:  a  coarser  grid  with  25  x  13  points  in  space,  40 
points  in  the  energy  to,  4  points  in  each  of  the  angular  /i  and  <f>  directions  (grid  B),  and  a 
finer  grid  with  91  x  33  points  in  space,  80  points  in  the  energy  to  and  12  points  in  each  of 
the  angular  /i  and  <fi  directions  (grid  C).  Figures  5  and  6  show  in  a  compact  way  the  scaled 
relative  differences  for  density,  electric  current,  electric  field  and  potential,  and  the  mean 
energy  multiplied  by  the  charge  density.  Measures  of  the  discrepancies  are  given  by  the 
following  formulas.  Let  q\y  be  the  value  of  an  hydrodynamical  quantity  qc  at  the  grid  point 
(. Xi,yj )  obtained  by  using  the  coarse  grid;  q{j  is  the  corresponding  value  obtained  by  using 
the  hire  grid.  Then, 


Qtj  -  4j 


dij  =  - jr-— -  (i  =  0, 

ma x{qj  }  —  nnnjg-'  } 

1, . v,,) 

and  (j  =  0,1  ,...,Ny) 

(2.14) 

A=  *  Yd,, 

l+N>h 

(*  =  0,1, 

(2.15) 

These  seem  to  be  good  functions  to  measure  the  dimensionless  distance  between  the  two 
surfaces.  Here,  Nx  +  1  and  Ny  +  1  are  the  numbers  of  the  fine  grid  points  in  the  x  and  y 
directions,  respectively.  We  use  a  simple  linear  interpolation  in  order  to  obtain  the  values 
of  the  quantity  qc  at  those  points  belonging  to  the  fine  grid  but  not  to  the  coarse  grid. 
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We  can  see  that  the  results  of  the  grid  A  and  grid  B  already  agree  reasonably  well,  even 
though  the  grid  B  has  only  as  few  as  4  points  in  some  directions.  The  relative  error  is  in 
general  less  than  5%,  except  for  the  region  near  y  =  0.2 ym  where  the  singularities  exist 
and  where  the  error  could  reach  around  10%  for  the  worst  resolved  quantity  (horizontal 
momentum).  This  indicates  that  the  code  can  be  used  on  very  coarse  meshes  with  fast 
turn-around  time  to  obtain  engineering  accuracy  with  average  relative  errors  below  5%. 
The  results  of  the  grid  A  and  grid  C  agree  remarkably  well  overall.  The  relative  error  is 
in  general  well  below  2%,  except  for  the  region  near  y  =  0.2 ym  where  the  singularities 
exist  and  the  error  could  reach  around  4.5%  for  the  worst  resolved  quantity  (horizontal 
momentum).  The  larger  errors  near  y  =  0.2 ym  are  caused  by  the  fact  that  the  singularities 
are  not  as  well  captured  with  the  coarser  grid  as  the  smooth  part  of  the  solution.  This 
can  be  seen  in  Figure  7.  This  is  particularly  important  from  a  computational  time  point 
of  view  since  the  fine  grid  results  have  taken  almost  2  weeks  to  achieve  steady  state  at  5ps 
while  the  coarse  grid  takes  less  than  an  hour  in  a  standard  Pentium  IV-2.4  GhZ  PC.  This 
extremely  high  computational  cost  of  the  fine  grid  is  due  to  the  restrictive  CFL  condition 
caused  by  singularities  more  than  the  increased  number  of  grid  points. 

In  order  to  reduce  the  computational  time  for  the  fine  mesh,  we  introduce  a  limiter  for 
the  electric  held  through  the  formula 

E Cut  =  (min{max{£'x?  —  Ec},  Ec}  ,  min{max{i%,  -Ec},  Ec}  ,  0)  (2.16) 

where  Ec  is  positive  and  fixed.  In  our  simulation  (grid  Ecut),  we  choose  Ec  =  150  (kV/cm) 
and  the  same  grid  of  the  case  A.  At  each  time  step,  we  solve  the  Poisson  equation  and 
evaluate  the  true  electric  held  E,  but  use  instead  Ecut  in  the  force  term  of  the  Boltzmann 
equation.  The  discrepancies  between  the  results  of  grid  A  and  Ecut  are  reported  in  Figure 
8  and  comparison  results  in  Figure  9  show  a  remarkable  agreement  with  respect  to  the 
original  results  at  a  much  lower  computational  time. 

2.5  Deterministic  simulation  of  a  2D  Si  MOSFET 

In  this  subsection  we  report  the  results  of  the  simulation  of  a  MOSFET  whose  geometry 
is  sketched  in  Figure  10  in  which  the  contacts  are  doped  with  an  electron  density  given 
by  n+  —  3  x  1017  cm-3.  No  doping  is  given  in  the  n  domain.  Here,  we  have  solved  the 
Boltzmann-Poisson  system  disregarding  completely  the  dynamics  of  the  holes.  The  origin 
of  the  potential  is  set  at  the  source  and  we  consider  applied  biases  of  Vgate  =  0.4V  at  the 
gate  and  Vdrain  =  IE  at  the  drain  respectively. 
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The  steady  state  results  are  reported  in  Figures  11  and  12  which  are  quite  clear  about 
the  qualitative  transport  of  electrons  through  the  device.  The  main  difference  with  respect 
to  the  MESFET  device  relies  on  the  fact  that  the  singularities  of  the  field  and  the  square- 
root  behavior  of  changing-data  points  on  the  electric  potential  are  in  the  oxide  region 
where  the  transport  of  electrons  is  not  computed  through  the  Boltzmann  equation.  Thus, 
the  appearance  of  the  singularities  of  the  force  field  does  not  affect  the  CFL  condition  in 
the  area  in  which  the  transport  equation  is  solved.  This  makes  MOSFET  much  easier  to 
compute  than  MESFET  in  terms  of  CPU  time.  The  results  shown  have  been  obtained  with 
a  grid  of  49  x  33  points  in  space,  66  points  in  the  energy  u  and  12  points  in  each  of  the  fi 
and  (j)  directions,  and  have  taken  3  days  on  a  standard  Pentium  1V-2.4  GhZ  PC. 

Similar  comments  with  respect  to  a  comparison  with  coarser  grids  can  be  made  in  the 
case  of  a  MOSFET.  Although  we  do  not  show  any  figure  to  save  space,  coarse  to  fine 
grid  comparisons  have  been  performed  and  the  conclusions  of  the  MESFET  remain  valid. 
Thus,  coarse  grid  BTE  results  are  competitive  in  accuracy  versus  computational  time  with 
respect  to  DSMC  results  for  a  MOSFET  device.  Finally,  we  report  in  Figure  13  the  pdf’s 
<F,  averaged  over  q f>,  at  different  points  of  the  device  showing  again  the  overall  smoothness 
of  the  distribution  function. 


3  Comparisons  with  Monte  Carlo  results 

The  classical  widespread  method  to  investigate  transport  processes  in  semiconductor  devices 
at  the  microscopic  level  is  based  on  the  Direct  Simulation  Monte  Carlo  (DSMC)  method. 
The  main  advantage  of  this  approach  consists  in  the  simple  way  to  take  into  account  different 
forces  acting  on  the  electrons.  The  dynamics  of  particles  is  well  described  by  Newton’s  law 
when  they  move  freely,  driven  by  the  mean  electric  field,  and  by  a  stochastic  process  during 
the  collisions  with  phonons. 

Although  DSMC  is  a  stochastic  numerical  method  and  the  BTE  equation  is  a  determin¬ 
istic  equation,  it  has  been  shown  in  [26]  that,  for  the  classical  case  of  space  homogeneous 
mono-atomic  elastic  gases,  the  stochastic  numerical  method  converges,  in  a  probabilistic 
sense,  to  the  solution  of  the  BTE.  This  fact  has  led  to  a  common  acceptance  on  the  equiv¬ 
alence  of  stochastic  approaches  to  solve  the  BTE  and  their  corresponding  deterministic 
approximations,  even  in  cases  of  strong  space  inhomogeneities. 

Recent  comparisons  for  Boltzmann-Poisson  transport  models  [18]  show  an  excellent 
agreement  when  spatially-homogeneous  solutions  (bulk  case)  are  analyzed  and  its  relevant 
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moments  are  compared.  This  confirms  that  collision  and  drift  (for  constant  electric  field) 
are  well-described  by  both  models.  We  stress  that  the  same  physical  parameters  were  used 
in  both  models  and,  in  particular,  the  scattering  rates  coincide  exactly.  In  particular,  when 
inhomogeneous  solutions  [3,  18]  are  considered,  the  agreement  between  hydrodynamical 
quantities  is  usually  good  for  small  relative  energies,  but  not  always  excellent. 

In  order  to  understand  the  origin  of  these  discrepancies,  it  is  useful  to  recall  two  aspects 
of  DSMC  methods.  The  first  one,  which  relates  to  the  boundary  conditions,  is  the  rule 
that  is  used  when  an  electron  must  be  injected  into  the  device.  The  second  one  relates  to 
the  length  of  the  free  flight  time.  Since  there  is  not  a  unique  resolution  to  these  two  issues 
[15,  13]  different  choices  give  often  inconsistent  results  when  comparisons  that  involve  these 
type  of  stochastic  methods  are  made.  These  facts  are  well-known  in  the  DSMC  community. 
However,  we  find  it  useful  to  discuss  these  two  factors  in  more  detail  in  order  to  make 
clear  the  comparison  criteria  between  stochastic  DSMC  and  deterministic  BTE  numerical 
solutions  and  to  explain  the  source  of  their  differences. 

We  first  consider  the  boundary  issue  for  the  stochastic  DSMC  method.  Since  the  elec¬ 
trons  can  exit  from  the  interior  of  the  device  (for  instance,  through  the  drain  or  source 
contacts  in  the  MESFET  or  MOSFET  geometry),  the  total  density  of  charges  may  de¬ 
crease.  Hence,  we  must  decide  whether  to  inject  electrons  to  balance  this  loss.  For  spatially 
homogeneous  solutions  where  charge  must  be  conserved,  this  is  a  simple  problem  to  ad¬ 
dress.  However,  for  general  geometries  with  strong  space  inhomogeneities,  the  difficulty  in 
treating  the  electron  injection  issue  arises  since  electrons  flow  through  metallic  contacts, 
which  are  also  connected  to  a  circuit.  In  principle,  a  full  treatment  of  the  system  consisting 
of  both  the  device  and  the  circuit  would  solve  the  problem  completely.  Such  simulation, 
however,  can  be  attained  only  in  very  simple  cases. 

Therefore  we  must  resort  to  some  reasonable  rule  which  might  be  justified  by  means 
of  asymptotic  analysis,  physical  considerations  or  empirical  tests.  In  this  context,  a  well- 
accepted  starting  point  to  decide  how  many  particles  must  be  injected  is  called  charge 
neutrality.  In  practice,  after  every  time  step  of  DSMC,  we  consider  the  cells  in  the  vicinity 
of  the  contacts,  usually  those  including  boundaries,  where  it  is  possible  to  inject  electrons, 
and  we  evaluate  the  charge  density  in  each  of  these  cells. 

In  particular,  we  End  that  the  charge  density  may  be  equal,  greater  or  smaller  than 
the  doping  density  on  every  cell.  In  the  first  case  no  change  is  required;  in  the  second  we 
easily  eliminate  a  few  particles  in  order  to  get  the  right  balance;  in  the  last  case,  we  have  to 
supply  electrons  to  balance  the  loss  of  charge.  This  electron  injection  represents  a  crucial 
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point,  since  we  must  decide  how  to  choose  the  velocity  of  the  entering  particles  (position 
can  be  chosen  randomly  inside  every  cell  due  to  the  stochastic  nature  of  the  simulation.) 

There  is  no  unique  recipe  for  electron  injection  in  order  to  balance  charges.  Simple 
common  chosen  rules  [13]  are  based  on  the  following  distribution  functions  or  methods: 

a)  The  global  Maxwellian  distribution  exp(— e/ksTi). 

b)  A  shifted  semi-Maxwellian. 

c)  A  velocity-weighted  Maxwellian  distribution. 

d)  The  reservoir  method. 

For  the  first  choice,  the  energy  of  the  entering  particles  is  given  according  to  the 
Maxwellian  distribution  function  and  then  the  magnitude  of  the  wave- vector  |k|  is  de¬ 
termined.  So,  in  order  to  fix  k,  we  need  the  values  of  the  two  angular  variables  which  can 
be  randomly  chosen  with  the  only  constraint  that  the  velocity  of  the  electron  is  directed 
towards  the  interior  of  the  device.  This  rule  is  reasonable  when  the  mean  current  is  almost 
zero  near  the  contact  which  is  the  case  for  highly  doped  regions  where  the  Debye  length 
asymptotics  dominates  the  transport  due  to  field  effects. 

Concerning  the  case  b),  we  must  point  out  that  the  meaning  of  the  shifted  Maxwellian 
in  the  framework  of  semiconductor  transport  differs  from  those  of  classical  perfect  gases. 
If  we  assume  the  parabolic  band  approximation  for  the  energy  band  structure,  then  the 
wave-vector  k  can  be  seen  as  a  classical  particle  velocity  and  a  term  as  k  —  u  (u  being 
a  fixed  mean  velocity)  is  a  shifted  velocity.  Moreover,  we  can  define  also  the  electron 
temperature.  However,  caution  is  needed  since  a  shifted  Maxwellian  is  not  a  solution  of  the 
BTE  except  for  the  case  in  which  only  electron-electron  interactions  are  considered.  This 
means  electron-phonon  interactions  are  neglected  and  that  the  force  field  term  in  the  free 
streaming  operator  takes  into  account  only  electron- doping  interaction  and  the  external 
electric  field.  When  a  non-parabolic  band  approximation  is  employed,  we  may  continue 
the  use  of  a  shifted  velocity,  but  the  physical  meaning  of  the  factor  inside  the  exponential 
function  is  not  clear  since  it  does  not  represent  the  kinetic  temperature  anymore.  To 
overcome  this  problem,  often  the  lattice  temperature  is  used  under  the  assumption  that 
contacts  are  in  thermal  equilibrium  with  the  crystal.  Therefore,  choice  b)  allows  us  to 
have  a  particle  moving  according  to  the  mean  electric  current,  but  with  an  inaccurate 
energy.  In  fact,  a  shifted  Maxwellian  when  used  in  hot  electron  situations  such  as  strong 
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inelastic  scattering  overestimates  significantly  the  energy  loss  from  electrons  to  phonons  and 
therefore  underestimates  the  electron  mean  energy.  Thus,  the  use  of  this  assumption  on  the 
incoming  distribution  is  reasonable  only  when  electron-electron  scattering  is  the  dominant 
process.  This  is  the  case  when  the  electron  density  is  large  enough  (a  value  greater  than 
5  •  101'  cm-3  seems  reasonable)  and  the  electric  held  is  not  very  strong. 

Some  comparisons  between  the  BTE  and  DSMC  results  exhibit  this  effect  in  ID  diode 
simulations  of  channels  under  400nm  subject  to  relatively  strong  voltage  bias.  In  reference 
[3]  we  implemented  choice  a).  The  agreement  is  excellent  except  in  the  case  of  a  50- 
nm  channel  diode,  where  a  small  difference  in  the  curves  concerning  the  hydrodynamical 
energy  near  the  extreme  of  the  diode  is  evident.  In  [18]  the  stochastic  Damocles  code 
was  employed  with  a  shifted  Maxwellian.  The  DSMC  and  deterministic  simulation  results 
agree  well  for  longer  devices,  but  as  the  length  of  the  channel  is  shorter,  the  behavior  of 
the  hydrodynamical  variables  differ  clearly. 

The  case  c)  employs  a  distribution  function,  which,  apart  a  normalizing  factor,  is  the 
product  of  the  global  Maxwellian  (with  the  lattice  temperature)  by  the  electron  velocity. 
In  ref.  [13]  an  interesting  application  to  a  ID  Schottky-barrier  diode  shows  a  good  perfor¬ 
mance  of  this  choice.  We  note  that  all  the  cases  a)-c)  assume  that  the  temperature  at  the 
contacts  is  very  near  to  the  lattice  temperature  and  the  shape  of  the  distribution  function 
of  the  entering  particles  is  given.  In  our  2D  simulations  the  appearance  of  the  singularities 
near  the  extremes  of  the  contacts  suggests  that  the  hypotheses  concerning  the  temperature 
might  not  be  accurate  at  these  points.  Hence,  we  believe  that  distributions  a-c)  are  not 
always  reasonable  in  these  cases.  Consequently,  we  employed  the  rule  d)in  this  paper.  We 
add  a  layer  of  ghost  cells  outside  the  device  near  the  drain  and  source  contact  regions.  Elec¬ 
trons  move  inside  these  cells  under  the  electric  held  obtained  by  the  simplest  extrapolation 
formula  using  the  nearest  values  of  the  true  cells.  Moreover,  electrons  encounter  collisions 
with  phonons. 

We  then  apply  the  charge  neutrality  criterion  only  in  the  ghost  cells  and  create  new 
particles  if  necessary  according  to  the  rule  a).  In  this  way  electrons  moving  from  the  ghost 
cells  to  the  device  enter  with  reasonable  velocity  and  energy.  In  our  DSMC  simulation,  this 
mechanism  seems  to  perform  better  than  the  other  approaches,  even  though  the  electric 
current  near  the  contacts  is  not  well  described  (see  Figures  15  and  19).  One  of  the  reasons 
for  this  failure  arises  from  the  number  of  ghost  cells  added  to  the  physical  device.  In  fact,  a 
small  number  guarantees  the  charge  neutrality  rule  at  the  contacts,  but  the  electric  current 
might  suffer  from  the  low-energy  entering  particles.  On  the  other  hand,  by  increasing  the 
number  of  the  ghost  cells,  we  improve  the  accuracy  of  the  electric  current  but  we  get  a 
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nonphysical  behavior  of  the  density  profile.  Finally,  we  also  prefer  to  use  rule  d)  because  it 
is  the  most  consistent  with  the  boundary  conditions  which  we  use  in  the  BTE  simulations. 

The  accuracy  of  the  deterministic  simulations  may  be  a  useful  tool  to  check  the  validity 
of  the  four  rules  a-d)  in  the  case  of  2D  devices.  Of  course,  this  study  requires  a  large  number 
of  simulations,  using  different  values  for  the  applied  voltages,  doping  and  geometry  sizes  in 
order  to  reach  meaningful  conclusions.  This  interesting  research  is  out  of  the  scope  of  the 
present  work. 

The  second  item  mentioned  above,  i.e.,  the  free  flight  duration,  is  based  on  the  total 
scattering  rate.  There  is  a  theoretical  formula  to  evaluate  this  parameter. 

We  denote  by  P[k(t)]dt  the  probability  that  an  electron  in  the  state  k  encounters  a 
collision  during  the  time  interval  dt.  When  considering  the  motion  of  this  single  electron, 
k(t)  denotes  the  wave-vector  k  as  a  function  of  the  time  t.  Let  t  =  0  be  the  time  of  the 
last  collision.  Then,  the  probability  V(t)  that  the  electron  will  encounter  its  next  collision 
during  dt  is  given  by 

V(t)dt  =  P[k(t)]  exp  —  f  P[k(t')]dt'  dt .  (3.1) 

L  Jo 

We  remark  that  we  must  use  the  above  formula  after  every  collision  and  for  any  particle, 
that  is,  if  we  consider  another  electron,  then  its  motion  will  be  different  and  therefore  it  will 
have  a  different  probability  P[k(f)]df.  Moreover,  it  is  evident  that  this  procedure  requires 
the  knowledge  of  the  function  P[k(f)]  as  t  — >  k(f).  This  whole  procedure  is  too  costly  for 
practical  applications,  and  thus,  simple  approximations  of  the  above  formula  have  been 
proposed  (see  for  instance  [15]).  In  our  simulations,  we  use  the  simplest  rule  given  in  [25]: 
fix  the  maximum  value  of  the  electron  energy  emax,  and  define 

T  =  max  |  [  S^k,  k')  dk'  :  e(k)  <  emax 

k  l  Jr3 

Then  approximate  P[k(t)]  from  (3.1)  with  the  constant  T.  So,  we  obtain 

V  =  Texp(-tT) . 

Since  r  =  P/T  is  a  random  number  distributed  uniformly  between  0  and  1,  we  can  choose 
the  free  flight  duration  tr  according  to  the  rule 

U  =  ~  y  log(r) . 

Therefore,  every  electron  has  a  proper  time  tr  which  is  again  generated  after  every  collision 
drawing  the  random  number  r. 
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Numerical  solutions  obtained  by  using  Boltzmann-Poisson  system  and  DSMC  are  con¬ 
sidered  in  the  stationary  regime  where  noise  of  the  stochastic  solutions  are  reduced.  In 
Figure  14,  the  error  formula  (2.14)  is  employed  to  compare  results  between  the  BTE  (grid 
A)  and  DSMC  simulations  for  MESFET  device.  We  also  report  BTE-DSMC  comparison 
results  of  three  cuts  for  fixed  values  of  y  —  0.2 fim,  y  =  0.175 ym  and  y  =  0.1  ym  in  Figures 
15,  16  and  17  respectively. 

These  results  show  excellent  agreement  between  the  deterministic  and  the  DSMC  simu¬ 
lations,  except  in  almost-empty  regions  (see  Figure  15).  Here,  we  are  only  stressing  different 
values  of  the  energy  near  the  gate-boundary.  The  mean  energy  for  DSMC  is  lower  than  the 
corresponding  value  obtained  by  the  deterministic  result.  We  believe  that  this  fact  depends 
on  the  energy  of  the  injected  electrons. 

Moreover,  energy  curves  corresponding  to  the  DSMC  simulation  show  unusual  non  phys¬ 
ical  oscillations  at  y  =  0.1  ym  near  the  boundary  x  =  0.  Here  the  electric  held  is  low  and  we 
believe  that  the  deterministic  result  is  very  accurate.  This  seems  to  indicate  that  in  these 
regions  the  free  flight  duration  is  inaccurate;  hence  electrons  do  not  have  a  correct  number 
of  collisions  with  phonons,  and  the  amount  of  their  energy  transferred  to  the  lattice  is  also 
inaccurate.  Analogous  conclusions  hold  also  in  the  case  of  a  MOSFET;  see  Figures  18,  19, 
20  and  21. 

Both  for  the  MESFET  and  the  MOSFET  device,  the  discrepancies  between  the  DSMC 
and  the  deterministic  results  are  important  near  the  ohmic  contacts.  It  is  reasonable  that  a 
more  refined  Monte  Carlo  code  is  able  to  generate  better  hydrodynamical  profiles.  In  this 
paper,  we  have  focused  our  attention  on  solving  the  BTE  equation  deterministically  and 
on  using  the  DSMC  scheme,  based  on  the  code  described  in  [25],  to  check  the  accuracy, 
the  performance  of  the  deterministic  numerical  scheme  and  the  agreement,  except  for  small 
spatial  regions,  with  stochastic  solutions.  We  remark  that  the  excellent  agreement  for  grid 
points  where  the  solution  shows  strong  gradients  confirms  the  absence  of  viscosity  effects 
in  WENO-scheme  results. 


4  Concluding  remarks 

We  have  adapted  a  high  order  accurate  finite  difference  weighted  essentially  non-oscillatory 
(WENO)  solver  for  the  direct  numerical  simulation  of  two  dimensional  Boltzmann  transport 
equation  (BTE)  coupled  with  a  Poisson  equation  describing  two  dimensional  semiconduc¬ 
tor  devices,  including  the  MESFET  and  MOSFET  devices.  We  have  demonstrated  the 
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performance  of  this  solver  through  a  comparison  between  coarse  and  fine  grid  simulations 
and  through  a  comparison  with  the  results  obtained  by  a  direct  simulation  Monte  Carlo 
(DSMC)  solver.  We  have  clarified  several  issues  related  to  the  two  dimensional  devices  in¬ 
cluding  the  implementation  of  boundary  conditions,  singularities,  and  factors  important  for 
the  comparison  between  deterministic  BTE  and  stochastic  DSMC  solvers.  The  simulation 
results  demonstrate  the  superior  capability  of  the  WENO-BTE  solver,  which  can  be  used 
to  obtain  accurate  results  with  a  very  coarse  mesh  (with  as  few  as  4  grid  points  in  the 
angular  directions). 
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Figure  2:  MESFET  device.  BTE  results  by  grid  A.  t  =  5  ps.  Top  left:  the  charge  density  p 
(10-17  •  cm-3);  top  right:  the  energy  £  (eV);  bottom  left:  the  electric  potential  V;  bottom 


right:  |Ej  ( kV/cm ). 
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Figure  3:  MESFET  device.  BTE  results  by  grid  A.  t  =  5  ps.  Left:  the  current;  right:  the 
electric  held  E. 


steady  state  pdf  in  (0.185,0.2) 


steady  state  pdf  in  (0.414,0.2) 


Figure  4:  MESFET  device.  BTE  results,  t  —  5  ps.  Pdf’s  in  steady  state  averaged  over 
(j)  at  different  locations  of  the  MESFET.  Top  left:  at  (x,y)  =  (0.185,0.2);  top  right:  at 
(x,y)  =  (0.414,0.2);  bottom  left:  at  (x,y)  =  (0.3,  0.1);  bottom  right:  at  (x,y)  =  (0.3,  0.2). 
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Figure  5:  MESFET  device.  BTE  results,  t  —  5  ps.  Discrepancies  (2.14)  between  the 
grids  A-B  as  a  function  of  y.  Left:  the  density  p  (continuous  line)  and  the  total  energy 
p£  (dashed  line);  right:  the  x- momentum  pux  (continuous  line)  and  the  ^/-momentum  puy 
(dashed  line). 
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Figure  6:  MESFET  device.  BTE  results,  t  —  5  ps.  Discrepancies  (2.14)  between  the 
grids  A-C  as  a  function  of  y.  Left:  the  density  p  (continuous  line)  and  the  total  energy 
p£  (dashed  line);  right:  the  x- momentum  pux  (continuous  line)  and  the  ^/-momentum  puy 
(dashed  line). 


32 


density  at  y  =  0.200  x-component  of  the  electric  field  at  y  =  0.200 


y-component  of  the  electric  field  at  y  =  0.200 


y-component  of  the  momentum  at  y  =  0.200 


Figure  7:  MESFET  device.  BTE  results,  t  —  5  ps.  Comparisons  between  different  grids  at 
y  =  0.2.  Solid  line:  grid  A,  dashed  line:  grid  C.  Top  left:  density  p;  top  right:  x-component 
of  the  electric  held  E;  bottom  left:  ^/-component  of  the  electric  held  E;  bottom  right:  the 
x-momentum  pux. 
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Figure  8:  MESFET  device.  BTE  results,  t  —  5  ps.  Discrepancies  (2.14)  between  the 
grids  A-Ecuj  as  a  function  of  y.  Left:  the  density  p  (continuous  line)  and  the  total  energy 
pS  (dashed  line);  right:  the  x- momentum  pux  (continuous  line)  and  the  ^-momentum  puy 
(dashed  line). 
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density  at  y  =  0.200  x-component  of  the  electric  field  at  y  =  0.200 


y-component  of  the  electric  field  at  y  =  0.200  y-component  of  the  momentum  at  y  =  0.200 


Figure  9:  MESFET  device.  BTE  results,  t  —  5  ps.  Comparisons  between  different  grids 
at  y  =  0.2.  Solid  line:  grid  A,  dashed  line:  grid  Ecu^.  Top  left:  density  p]  top  right:  x- 
component  of  the  electric  held  E;  bottom  left:  y-component  of  the  electric  held  E;  bottom 
right:  the  ^-momentum  pux. 
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Figure  11:  MOSFET  device.  BTE  results,  t  —  5  ps.  Top  left:  the  charge  density  p 
(10”1'  •  cm-3);  top  right:  the  energy  £  (eV)]  bottom  left:  the  electric  potential  V ;  bottom 
right:  |E|  ( kV/cm ). 
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Figure  12:  MOSFET  device.  BTE  results,  t  —  5  ps.  Left:  the  current;  right:  the  electric 
field  E. 
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steady  state  pdf  in  (0.04,0.2)  steady  state  pdf  in  (0.24,0.12) 


steady  state  pdf  in  (0.24,0.2) 


Figure  13:  MOSFET  device.  BTE  results,  t  —  5  ps.  Pdf’s  in  steady  state  averaged  over 
(f)  at  different  locations  of  the  MOSFET.  Top  left:  at  (x,y)  =  (0.04,0.2);  top  right:  at 
(x,y)  =  (0.24,0.12);  bottom:  at  (x,y)  =  (0.24,0.2). 
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Figure  14:  MESFET  device,  t  —  5  ps.  Discrepancies  (2.14)  between  the  BTE  result  (grid 
A)  and  the  DSMC  result  as  a  function  of  y.  Left:  the  density  p  (continuous  line)  and 
the  total  energy  pS  (dashed  line);  right:  the  ^-momentum  pux  (continuous  line)  and  the 
{/-momentum  puy  (dashed  line). 
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x-component  of  the  momentum  at  y  =  0.200  y-component  of  the  momentum  at  y  =  0.200 


Figure  15:  MESFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at 
t  —  5  ps  and  y  =  0.2.  Solid  line:  the  BTE  results  (grid  C),  dashed  line:  the  DSMC  results. 
Top  left:  density  p\  top  right:  energy  £:  bottom  left:  the  ^-momentum  pux\  bottom  right: 
the  y- momentum  puy. 
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density  at  y  =  0. 175  energy  at  y  =  0. 175 


x-component  of  the  momentum  at  y  =  0. 175 


y-component  of  the  momentum  at  y  =  0. 175 


Figure  16:  MESFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at  t  —  5 
ps  and  y  =  0.175.  Solid  line:  the  BTE  results  (grid  C),  dashed  line:  the  DSMC  results. 
Top  left:  density  p;  top  right:  energy  £;  bottom  left:  the  ^-momentum  pux\  bottom  right: 
the  y- momentum  puy. 
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density  at  y  =  0. 100  energy  at  y  =  0. 100 


x-component  of  the  momentum  at  y  =  0. 100 


y-component  of  the  momentum  at  y  =  0.100 


Figure  17:  MESFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at 
t  —  5  ps  and  y  =  0.1.  Solid  line:  the  BTE  results  (grid  C),  dashed  line:  the  DSMC  results. 
Top  left:  density  p\  top  right:  energy  £:  bottom  left:  the  ^-momentum  pux\  bottom  right: 
the  y- momentum  puy. 
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Figure  18:  MOSFET  device,  t  —  5  ps.  Discrepancies  (2.14)  between  the  BTE  and  the 
DSMC  results  as  a  function  of  y.  Left:  the  density  p  (continuous  line)  and  the  total  energy 
p£  (dashed  line);  right:  the  x- momentum  pux  (continuous  line)  and  the  ^/-momentum  puy 
(dashed  line). 


44 


density  at  y  =  0.240  energy  at  y  =  0.240 


x-component  of  the  momentum  at  y  =  0.240  y-component  of  the  momentum  at  y  =  0.240 


Figure  19:  MOSFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at 
t  —  5  ps  and  y  =  0.24.  Solid  line:  the  BTE  results,  dashed  line:  the  DSMC  results.  Top 
left:  density  p\  top  right:  energy  £ ;  bottom  left:  the  x- momentum  pux\  bottom  right:  the 
y- momentum  puy. 
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density  at  y  =  0.200  energy  at  y  =  0.200 


x-component  of  the  momentum  at  y  =  0.200  y-component  of  the  momentum  at  y  =  0.200 


Figure  20:  MOSFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at 
t  —  5  ps  and  y  =  0.12.  Solid  line:  the  BTE  results,  dashed  line:  the  DSMC  results.  Top 
left:  density  p\  top  right:  energy  £ ;  bottom  left:  the  x- momentum  pux\  bottom  right:  the 
y- momentum  puy. 
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density  at  y  =  0. 1 20  energy  at  y  =  0. 1 20 


x-component  of  the  momentum  at  y  =  0.120 


y-component  of  the  momentum  at  y  =  0. 120 


Figure  21:  MOSFET  device.  Comparisons  between  the  BTE  and  the  DSMC  results  at 
t  —  5  ps  and  y  =  0.12.  Solid  line:  the  BTE  results,  dashed  line:  the  DSMC  results.  Top 
left:  density  p\  top  right:  energy  £ ;  bottom  left:  the  x- momentum  pux\  bottom  right:  the 
y- momentum  puy. 
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